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ABSTRACT 

■^ , A kinetic equation for the collisional evolution of stable, bound, self gravitating and 

r •^ ' slowly relaxing systems is established, which is valid when the number of constituents 

is very large. It accounts for the detailed dynamics and self consistent dressing by 
collective gravitational interaction of the colliding particles, for the system's inhomo- 
geneity and for different constituent's masses. It describes the coupled evolution of 
O ' collisionally interacting populations, such as stars in a thick disk and the molecular 

4— > , clouds off which they scatter. 

^ ' The kinetic equation derives from the BBGKY hierarchy in the limit of weak, but 

non- vanishing, binary correlations, an approximation which is well justified for large 
stellar systems. The evolution of the one-body distribution function is described in 

" t-j^ , action angle space. The collective response is calculated using a biorthogonal basis of 

1/-^ ' pairs of density-potential functions. 

Q^ , The collision operators are expressed in terms of the collective response function al- 

lowed by the existing distribution functions at any given time and involve particles in 
resonant motion. These equations are shown to satisfy an H-theorem. Because of the 

lO ' inhomogeneous character of the system, the relaxation causes the potential as well as 

^^ , the orbits of the particles to secularly evolve. The changing orbits also cause the angle 

Fourier coefficients of the basis potentials to change with time. We derive the set of 
equations which describes this coupled evolution of distribution functions, potential 
and basis Fourier coefficients for spherically symmetric systems. In the homogeneous 

K> , limit, which sacrifices the description of the evolution of the spatial structure of the 

\^ ' system but retains the effect of collective gravitational dressing, the kinetic equation 

reduces to a form similar to the Balescu-Lenard equation of plasma physics. 

Key words: stellar dynamics-galaxies: star clusters-plasmas 



1 INTRODUCTION AND MOTIVATION 

The description of collisional relaxation in a self-gravitating system usually rests on a Fokker-Planck equation in which the 
diffusion and braking coefficients are calculated in the local approximation, taking the finite dimensi on of the system into 



account by limiting the impact parameter of the collisions to a length of order of the system's size (jChandrasekhail 11942 



19431 : iBinnev fc Tremainelll987l : ISpitzerlll987l ). Although characteristic relaxa tion times may be somewhat overestimated by 



this approximation due to the neglect of collective self-gravitational effects ( Weinberg! Il993l ) . such a kinetic equation may 
provide in practice a reasonable description of the collisional relaxation of gravitationally bound systems. It nevertheless rests 
on assumptions which, from a principle point of view, are unsatisfactory because the motion of particles during the coUision is 
regarded as rectilinear and uniform and the system's inhomogeneity, which is basically the reason why collisions with an infinite 
impact parameter do not occur, is treated by way of an ill-defined cutoff. Moreover, the collective response of the system 
is not taken into account, since the Fokker-Planck collision term only considers binary collisions between naked particles. 
A self-gravitating medium, unlike an electrical plasma, does not respond to the presence in it of a particle by screening its 
interaction potential with other particles. As a result, even distant particles effectively interact, while in electrical, globally 
neutral, plasmas, the effective interaction distance is limited to the Debye length. In a self-gravitating system, the distance 
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between interacting particles is only limited by the system's inhomogeneity. The spatial structure of the system matters as 
well as the details of the particle orbits. 



The consistent inclusion of collective screening effects in a kinetic equation fo r electrically in tera cting weakly c oupled 
particles has been one of the major theoretical achievements in plasma physics w hen Balescul ( 1960l ) and iLenardI ( 1960f ) could 
derive an equation surpassing in consistency the simple Fokker-Planck equation ( Spitzeiill962l ). It is the aim of this paper to 
derive a similar equation for self-gravitating systems. The task is slightly more difficult because the screening of the electrical 
interaction at the, usually small, Debye length allows, in electrically interacting systems, to take the homogeneous and uniform 
motion limits. These limits cannot be taken in a self-gravitating system. We overcome this difficulty by expressing the kinetic 
equation in action angle space rather than in position momentum space. This is possible when the Hamiltonian corresponding 
to the average potential U{r) of the system is integrable. It is nevertheless uneasy in general to toggle from one to the 
other space, although this is certainly possible for spherically symmetric potentials, for flat systems (which may however be 
unstable) and for special thick disk potentials. Nu merical methods could be used to achieve the necessary transformation 
( Pichon fc CannorJll997l : [McMillan fc Binnevll2008h . As an illustrative example, we shall give special attention to spherically 
symmetric potentials, expanding their kinetic equation into a system which almost entirely avoids any calculation in the 
position-momentum space. The system's inhomogeneity requires that solutions to the Poisson equation are easily found for 
any inhomogeneous mass distributions. This is achieved by projecting on a biorthogonal basis of pairs of density-potential 
functions. 



Many astrophysical systems which have evolved to a quasi-stationary coUisionless equilibrium still keep evolving on time 
scales longer than the dynamical time as a result of gravitational noise induced by their own co nstituents or by e xternal ones. 
We disregard external perturbators, which we define as unbound to the system, although, as did lWeinbergj ( 2001a ). these could 
be treated, if numerous and frequent enough, as a given, non-evolving, population providing a source of gravitational noise for 
other populations. Loosely bound satellites or remote star populations are regarded as internal to the system. This is possible 
because our set of kinetic equations allows to simultaneously follow different mass populations. Dwarf satellite galaxies could 
be regarded for example as one such mass population. Globular clusters, dwarf galaxies, disk galaxies and their haloes are 
examples of bound systems still evolving as a result of internal noise caused by particle discreteness. Such systems are the 
object of our study. As in any weakly coupled system, the particles suffering collisions are dressed by the polarization clouds 
caused by their own infiuence on o ther particles. C ollisions between dressed particles have quantitatively different outcomes 
than collisions between naked ones ( Weinberdll998l ). This may reflect in significant differences in calculated effective relaxation 
times and braking or diffusion coefficients, especially when the system, though stable, is not too far from instability ( Weinberd 
1993) . It is therefore useful to account for collective dressing when calculating such processes as secular thick disk evolution, 
mass segregation in galaxies or in star clusters, or the damping by dynamical friction of galactic populations on high energy 
orbits. For simplicity, the kinetic equations to be derived below assume that the system is stationary on a dynamical time 
scale. They thus cannot address questions in which the distribution in angle variable matters, such as the dissolution of 
freshly accreted satellites, although a simple extension of the theory could. Since however our equations describe the coupled 
evolution of all populations present in the system, they are well suited to study, for example, the simultaneous evolution by 
dynamical friction and diffusion of a stellar population and the population of molecular clouds off which these stars scatter. 



The collective resp onse of a self gravitating system to the presence of a perturbing bo dy has been cons idered by a number 
of authors, analytically (IWeinberelll989l . ll995l : lMurah fc Tremainell 19981 : ISaha fc Jo3l2006l ) or numerically (JThielheim fc Wolfl 
1^84; Gncdi n fc Ostrikerll 19991 ). Sometimes, the reaction of this perturbation on the perturbing body itself is calculated, as did 
Kalnaia (,1972r i. who computed the drag on a large body mov ing in an homogeneous medium, taking the collective response of 
this medium into account, and lTremaine fc Weinberd ( 1984 1. who considered the global, self-consistent, perturbation caused 
by a satellite or a barred structure in a spherically symmetric system and its reaction on the perturbator object by the 
effect of dynamic al friction. The secular evolution of the system in response to su ch perturbations has bee n considered by 
Weinberd ( 2001al ). who considered general types of perturbations on a galaxy, and bv lPichon fc AubertI ( 20061 ) who considered 
perturbations caused by the cosmological environment on dark matter haloes. This evolution is of course in principle ob servable 



in N- body simulations, which however have their own difficulties in calculating the lo ng term evolution of such systems (JBinnev 



2004 ). A number of authors ( Muralilll999l : IWeinberdl2001al : |Pichon fc Auberdl2006l ) have studied the collective perturbati 
caused in a massive spherical galactic halo by its environment. They could calculate the respo nse of thi s syste m by resorting to 
a representation of the particle's motion in action and angle variables, a method first used bv lKalnaid (11977). We f ollow them 
on this road. They also made good use of a basis of biorthogonal pairs of density-potential functions. IWeinberd ( 19931 ) first 
derived a kinetic equation for the coUisional relaxation of a self gravitating system along these lines. His equation accounts 
for the self-consistent gravitational dressing of the particles, but is otherwise simplified, the geometry supposedly being that 
of an homogeneously filled periodic cube. The inhomogeneous nature of the system should be d escribed more ac curately, still 
accounting for collective gravitational dressing effects. This is specifically the aim of this paper. IChavanid (|2007l ) presented a 
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similar approach to ours for one-dimensional systems, the constituents of which interact by a general long range force. In this 
paper we further elaborate in section |8] on the structural evolution of the inhomogeneous system and on the secular evolution 
of the orbits. 



2 CUTTING THE BBGKY HIERARCHY 

2.1 Reduction of the hierarchy to a kinetic equation 

The Liouville equation for the N-body distribution function of a system of interacting particles can be translate d into a 



hierearchy of equations, the BB GKY hierarchy, for the reduced 1-body, 2-body, 3-body etc .. distribution functions (|Balescu 



19631 : iBinnev fc Tremainelll987l ). The equation for the 1-body distribution function also involves the 2-body distribution, the 



equation for the 2-body distribution involves the 3-body distribution and so on. The kinetic equation being meant to be an 
autonomous equation for the 1-body distribution /i(r, p, t), its derivation necessarily involves some approximation allowing to 
cut this hierarchy. This is usually done at the level of the equation of evolution of the 2-body distribution function, reducing 
it to a relation between the 2-body and the 1-body distribution functions. The simplified equation for the 2-body distribution 
/2(ri,pi, r2, P2, t) is then solved in terms of the 1-body distribution /i(ri,pi,t) and the result, once introduced in the first 
equation of the hierarchy, provides the desired kinetic equation for /i . 

Plasma physics knows of two such successful approximations: rare and short range interactions, allowi ng to ignore 3-body 



coUis ional effects on the evolution of the 2-body distribution function, leading to the Boltzmann equation (jUhlenbeck fc Ford 



19631) and weakly coupled, collective, systems in which the 3-bod y correlations may be neglected and the 2-body correlations 



considered weak, leading to the Balescu-Lenard equation (iBalescuiil960 : Lenard 19601 ). The weak correlation approximation 
is valid when the number of particles in the efi'ective interaction sphere, the Debye sphere, is large. This approximation is also 
valid for self-gravitating systems with a large number A'' of simultaneously interacting particles. The coupling in this case is 
indeed weak, the ratio of the average interaction energy to the average kinetic energy scaling as N~^'^ . This provides a solid 
basis for the derivation of a kinetic equation. The larger A*', the more valid the approximation is. For systems with a very 
large number of bodies, the resulting kinetic equation is almost exact, but for the description of strong collisions. 

The constituents of the system are considered to be point-like objects of difi'erent masses, which we refer to as particles. 
They need not all be stars, but could be other entities as well, such as molecular clouds, bound clusters, a population of 
satellites or lumps of dark matter in the halo of a galaxy. The kinetic equations to be derived below are valid as long as most 
collisions are weak, which implies that the coUisional evolution time of any type of particles remains long compared to the 
dynamical time. We assume that the masses of the consituents come in a finite set. Each mass group is labeled by a lower 
case latin letter. 

2.2 Notations 

An efficient and concise notation is needed. Some weakly relevant variables, such as time, will often be omitted from the list 
of arguments of some functions. The subscripts 1 or 2 on one- or two-body distributions or correlation functions will also be 
omitted, the number of arguments indicating the number of bodies involved. The 1-body distribution function of particles 
of species a (that is, of mass rria) is denoted by /", the 2-body distribution function of a pair of particles of species a and 
h (where a and h may be equal or difi'erent) is /"'' and the corresponding 2-body correlation function is p"'' — f'' — f^ f^- 
The space and momentum integral of a 1-body distribution function is the total number of particles of the considered species. 
Similarly, the space and momentum integral of 2-body distribution functions is the total number of pairs of the considered 
species. When a = b, pairs should be regarded as ordered entities. The position and momentum (ri, pi) of a particle is simply 
noted 1, for brevity. The three angle and three action variables of this particle similarly form a pair of vectors (wi, Ji). The 
same shorthand notation, 1, is used where the context commands. The notation dl represents either cfiricfipi or cPwicP Ji. 
These phase space volume elements are equal because both sets of variables are canonical. The velocity of particle 1 is vi . 
The gradient with respect to a vectorial variable u, like r, p, w or J, is denoted by Vu. The derivative with respect to time 
is noted dt- 

G being Newton's constant, the gravitational force suffered by a particle of species a with dynamical variables 1 (that is 
at ri with momentum pi) from a particle of species h with dynamical variables 2 is: 

F„i,(l, 2) = Gm,m6 ""^''^.^ • (1) 

We ignore any external force, be it tidal or exerted by some closeby external body. The collective gravitational force FIJ(1) 
exerted at ri on a particle of species a is the 1-body and species average of Fa6(l, 2): 
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F°(l)=^ /'d2F„,(l,2)/(2) . (2) 

b •' 

The gravitational potential U{vi) from which this force derives is: 

C/(ri) = -V /d2^^/''(2) . (3) 

•^ — ' / r2 — ri 

2.3 Weak correlations in terms of one-body distributions 

The first equation of the BBGKY hierarchy can be written: 

9t,r(l)+vi.Vri,r(l) + F°(l).Vp,r(l) = -^ /'d2F,,(l,2).Vp,p'"'(l,2) . (4) 

Neglecting 3-body correlations, the second equation of the BBGKY hierarchy can be written as: 



a, 



:ff"''(l,2)+(vi.Vri+V2-Vr,)5'"'(l,2)+(F°(l).Vpi+F?(2).Vp,)<7''^(i,2) 



^ . d3 /^(2,3) F,,(l,3)-Vp,r(l) +2J / d3 g^^{l,i) F6c(2, 3)- Vp,/(2)= F,t,(l, 2) ■ (Vp, - VpJ/" (i)/(2)-(5) 

c c 

Equation ([5]) is linear in the correlation function and has on its right hand side a source term S°'^{1, 2, t) which is a functional 
of the 1-body distribution functions, namely: 

S'"'(l,2,t)=F„,(l,2).(Vp,-Vpjr(l)/''(2) . (6) 

The solution for (?"''(1, 2, t) can be found in terms of the sources S by working out the Green's function, or propagator, of the 
operator on the left hand side of equation ([5]). This Green's function is a matrix in particle species space, CJpg(l,2, l',2',r), 
in terms of which the correlation function can be expressed as: 

dTJdl'jd2' g;l{l,2,l',2\T)S^\l\2\t-T) ■ (7) 

p,q 

Equation ([7| expresses the correlation function as a functional g'^''{l, 2; /) of the 1-body distributions. Once the 2-body 
propagator has been found, the solution ^ for (?"''(1,2) may be substituted on the right hand side of equation Q, which 
then depends explicitly, and only, on the 1-body distributions. We call it the collision operator C"(/) for species a: 

C^if) = -Y. h ^'"'(^' 2) • Vpi<?"'(l, 2; /) • (8) 

b •' 

The initial value of the 2-body propagator is: 

g,t(l, 2, 1', 2', 0) = 5; 5l 5(1 - 1') 5(2 - 2') , (9) 

where 5(1 — 1') is a Dirac function and 5^ a Kronecker symbol. By substituting equation ((Tj) in equation ((5]), it can be shown 
that the 2-body propagator can be factored into the product of two 1-body propagators: 

ept(l,2,l',2',r) = ep"(l,l',r)e,'(2,2',r) • (10) 

Had we considered strong interactions as well, the correlation function (;"''( 1,2) would not have been negligible compared to 
/"(l)/''(2) and the right hand side term of equation ((5]) would have been changed by the substitution of /"Z'+g"'' to /"/''. 
The 2-body propagators would in this case not factor as in equation (|10p . In the weak correlation approximation considered 
here, the 1-body propagators Qp{\,\' ,t) satisfy the linearized Vlasov equations: 

a.gp"(l,l',r)-f (vi-Vri+F°(l)-Vpjep"(l,l',r)-f (^ /"d2 5^(2,1', r)F„,(l, 2) j ■ Vp,.r(l) = , (11) 

with initial condition Qp(l,l' ,Q) = S^Sil — 1'). The solution of equation llllll ha s to be found for r > only, because of 



causality. According to Bogoliubov's synchronisation hypothesis ( Bogoliubovlll946l l. the 1-body distribution functions can be 
regarded as constant in equations ([7| and (|lip because they evolve on the relaxation time scale, which is much longer than 
the time required for the correlation function to reach an equilibrium, given the present value of the 1-body distributions. 
The correlations at a given time t then are functionals of the one particle distribution functions at this very same time. 

Equation pif) can be solved by means of a Laplace transform with respect to the time lapse r. The Laplace transform 
fiui) of a function of time fit) depends on a complex argument uj. The transformation and its inverse are defined by: 
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f{u)^ j f{t)e'-'dt and f(t)^^jf{u^)e-^'^'du>- (12) 

The direct transform is convergent only when the imaginary part of uj exceeds some ordinate of convergence, above which 
the function /(oj) is regular. Below it, it is defined by analytical continuation. The Bromwich contour B which appears in 
the inverse transformation runs parallel to the real axis from — cxa to +00 above all singularities of f{ijj)- Equation (|11|) is 
Laplace-transformed into: 



.i^e;(i,i',^) + (vi.v., + vl{i)- Vpj g;{i,i' ,u) + Y^ fd2g;{2,i',uj) F,e(i,2)- Vp,r(i) = s;s{i~i') 



(13) 



3 PARTICLE MOTIONS AND BASIS FUNCTIONS IN ANGLE AND ACTION VARIABLES 

3.1 Angle and action variables 

The particle motions in the self gravitational field are complex in general. This precludes a direct solution of equation (lllf) 
by integration along unperturbed trajec tories. It is pre ferable to change the position and momentum variables for a set 
of canonical angle and action variables (Goldstein! '1956). So doing, the description of the motion becomes simple, all the 
complexity being embodied in the relation between position and momentum variables and angle and action variables. By 
definition, the Hamiltonian Ti. in angle and action variables depends only on the three actions Ji, J2, J3, which we regard as 
the three components of an action vector J. The three actions are constants of the motion and the three angles wi,W2,W3 
which similarly form the components of an angle vector w, vary linearly in time. The angular frequency of the angle Wi is 
rii — &H/dJi. The frequencies fii form the components of a frequency vector fl which depends on J. For brevity, we use 
shorthand notations, such as: 

ill = fl{Ji) , fl[ = Sl{J[) ■ (14) 

The derivative following the motion is (v ■ Vr + F" • Vp). The actions being first integrals, this operator translates in angle 
and action variables into {dew /dt) ■ Vw, that is, v • Vr + F° • Vp = fJ • Vw The last, collective, term of the left hand side of 
equation (|13l) must be expressed in action and angle variables. It is of the frequently met general form: 

A ld2M{2)--^^—^-\/p,N{l) ■ (15) 



The force in equation (|15|) can be expressed in terms of a "potential" (j> such that: 

d2M(2) , '''"''' ^-Vri0(ri) ■ (16) 

'r2 — ri l'^' 

This potential and the "mass distribution" D from which it derives, depend on the function M(2) only. They are defined by: 



0(ri) = - d2- ^, D{r2)= d% M{2) ■ (17) 

J I r2 - ri I J 

Since 0(1) is independent of pi, Vri';/>(1)- Vp^ A'^(l) is the Poisson bracket { (j!)(l), A^(l) }. This bracket being invariant on a 
change of canonical variables, the expression (|15|l can be written as: 

jd2M{2) ,^^~^\3 ■ VpiiV(l) = -Vri<^ ■ Vp,JV(l) = -(Vwi<^-VjiiV(l)-Vji<?i-VwiiV(l)) • (18) 

All functions depend periodically, with period 27r, on the angles, with respect to which a discrete Fourier transform can be 
made. All components of the associated wave vector k are relative integers. The transform of any function /(w, J) and the 
inverse transform are defined by: 

/(w,J) = ^/k(J)e'''"' and /^(J) =^ |i^ /(w, J) e"'"- • (19) 

Each integral in the second term of equation (|19p is over the 27r period of one of the components of w. The transform of the 
Dirac function (5(w) is l/Svr'^ and the transform of unity is (5(k), where 5 is here a triple Kronecker symbol. The position r of 
a particle is a function of its angle and action variables, w and J. The simple Fourier transforms with respect to the angles 
wi of i/)"(ri) and Gp(l, 1' ,ijj) and the double Fourier transform of the propagator with respect to angles wi and w'l are: 

r(l) ^ V'£i(Ji) e;(l,l',^) « Gl[{3,,l',uj) gp"(l,l',^) ^ G^^^,(Ji,J'i,c.) . (20) 
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3.2 Biorthogonal density-potential bases 

A basis of biorthogonal de nsity-potential pairs is effective in c a lculat i ng tlie potent i al <j)(l) defined by eg. (1171). Many such 
bases have been proposed (lKalnaislll97ll: Clutton-Brocklll972l. Il973l: lKalnaidll976l: Uoki fc Iy3ll97S : Uoki. Noguchi fc lye 



I979I : ISahalll99ll : iHernauist fc Ostrikeilll992l : iRobiin fc Earnlll996l : iBrown fc Papaloizoulll998l : iRhamati fc Jallalilboogf ). A 



basis element is labeled by a greek letter. The dummy index rule is used for these basis indices. Let D°'{r) and tp°'{r) be the 
density and the potential of the element a of the basis. The potential ip" derives from the density distribution D" and is 
related to it by: 



r(r) 



dV^^ 



|r'-r I 
The basis is biorthogonal and normalized, such that: 



d\ D"{r) (^/'''(r))* = - S: 



13 



(21) 



(22) 



The symbol on the right of equation (|22p is a generalized Kronecker. The minus sign results from the fact that when a = 13 
the left hand side of equation (|22p necessarily is negative. The functions to be expanded on the basis being real, the complex 
conjugates of D" and ^", 



D"(r) = (D°(r))* and V°(r) = (V'"(r))* 



(23) 



also form an element a of the basis, which in general is different from a. The variable r being a length and the Kronecker S 
in (|22|) being dimensionless, equations (|2H) and (|22p imply that D" and tp" have dimensions L~"''^ and L~^'^ respectively. 
Any density distribution D{r) and its associated potential 0(r) can be expanded on the basis as: 



I3(r) = a,D"(r) 



(r) = aatp^ir) 



(24) 



The basis functions ^"(r) are not real in general, which implies that i/'"]^ 7^ ('/'k)*- The notation ?/>£* denotes the complex 
conjugate of tp^. The notation i/i^ is adopted for the k-Fourier transform of the function tp"{r) = {tp°'{r)y. In general, 
'/'k 7^ V'k*- Complex conjugation implies however that: 



V'k = (i/'-k)* 



(25) 



The coefficients of the expansions (I24|l can be calculated by using the biorthogonality relation l|22|) and expressed in angle and 
action variables by using the density-potential basis and angle Fourier coefficients. In particular, the coefficient Ua associated 
with the density field of equation (|17|l is: 



- fd2 M(2)(V'"(2))* = -Stt^^ fd^'Ji Mk(J2) Vr(J2 



(26) 



The expression (|15p is transformed in angle and action variables by using the density-potential basis and angle Fourier 
coefficients into: 



A d2M{2) 



r2 - ri 

r2 — ri I' 



VpiiV(l) 



A > Oq e 



where the expansion coefficients aa are given in terms of the function A/(2) by equation 
established below for the case when M(2) is the propagator, as in equation p3[) . 



-™i(v£i(l)*i.Vj,7V(l)-(Vj,v^£,(l))-VwiiV(l)) , (27) 

Other expressions of aa will be 



4 THE LINEARIZED VLASOV PROPAGATOR 

In a relaxing system, a coUisionless equilibrium is supposedly reached on a time scale shorter than the relaxation time, 
so that the system is stationary on the dynamical time scale. This means that the distributions /°(1) really are functions 
/"(Ji) of the actions only. The third, "collective", term on the left of equation p3p is of the form displayed in equation 
(|15|l . The corresponding factor A and functions A'^(l) and M{2) particularize in this case to A = Graamc, N{1) — /"(l) and 
M{2) = Gp{2, l',u])- For these functions A''(l) and M(2), the coefficients Ua of equation PSJI are: 



a^^(l',a;) 



-87r'J2fd'j G^''(J,l',a;)Vr(J) 

1, -^ 



Species-cumulative coefficients A^ are defined by: 
AS(l',a;) = ^m,a^''(l',a;) • 



(28) 



(29) 
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Equation (|13|) for the 1-body propagators is Fourier transformed with respect to wi using equation (|27|) . which gives: 
G^^(J.,l',^) = ^5,"5(J,-j;)-£:^- Gm. %^^&^ <(Ji)^^.(l',-) • (30) 



The coefficients A (eq. (|29|) ) can be expressed in terms of the Fourier transform of the propagators by using equation (|28|) : 

AS(l',a;) = -87r'^^ /"d^J m,V£*(J)G;^''(J,l',^) • (31) 

c k "' 

Operating on equation H30p as on the function G in equation (|3ip . a linear system is obtained for the species-cumulative 
coefficients A. It can be written: 

e"''H^^(l',c^) = aS(l',a;), (32) 

'7q(1 ,^) = -^J^P > ^--pjT , (33 



£"''H = r''-^^ /dVi STT^Gm^ ^k:(l)<(l) 



VZiZil) . ,34, 

o; — ki -iii 



The solution of equation (|32|l . obtained by inverting the matrix £'^'^(0;), is then introduced in equation pO[) . giving the Fourier 
and Laplace transform of the propagator. So doing, a function T> appears in the solution, which is defined by: 

- — ] ^ ^£,(Ji)(£-^H)"%f;(j;) . (35) 

Performing the inverse Fourier and Laplace transforms of equation (|30|) . the 1-body propagator itself is eventually found: 

r-^n 1' ^ ( '^^-'-r ^^ ^e-CKi-w.-k^.w;) / 87r^Gm„mp(ki.Vj,.r(l)) \ 

^-^^'^^^^ = ]^2^' 2^2^ 8.3(,_k,.n0 (^'5p'5(ki-kO<5(Ji-Ji)+ (,_k;.o;)i,^^^,(j^,j,,^)j ' (36) 



5 THE KINETIC EQUATION 

5.1 Explicit writing of the kinetic equation 

The correlation function is obtained from the solution p6|l for the 1-body propagator by using equations pop and (O. The 
kinetic equation and its collision operator are then given by equations (Q and ([H]). Thanks to the Bogoliubov synchronisation 
hypothesis, this equation is local in time, because the source term S''"'(l', 2',i — r) in equation can be regarded as 
independent of r and equal to its value at t = 0. The collision operator for the evolution of the distribution function of species 
a, C'^{f), is defined by equation ((8]) and can be written as: 

••■ (F„,(l,2).VpJ {g;{l,l',u)g',{2,2',u') (F,,(1',2')- (Vp, -Vp,)) r(l',t)r(2',t)) . (37) 

The somewhat lengthy transformations that must be performed to express this equation in terms of the angle and action 
variables, of the density-potential basis and of its angle Fourier transforms are described in appendix]^ They eventually yield 
the following final form of the kinetic equations: 

ftr(Ji) = EEE^''^^ ^-'G'mlml k,.V., ( J^''""7'?'"!\,. (ki.V.,-k..V.,),r(Ji)/(J.)) , (38) 

bl^T^-^ VjOkik2(Ji,J2,ki-ni)| J 

where V is defined by equation H35p and the response matrix elements e"'^ needed to determine V are expressed in terms of 
the 1-body distribution functions by equation (|34p . No convective term f2i-Vwi/"(l) appears on the left hand side of eq.i 
because in a slowly relaxing system the distribution functions /"(I) are meant to depend only on the actions. 



5.2 Physical content of the kinetic equation 

Equation (|38p describes the relaxation of the distribution functions caused by the, supposedl y weak, noise created by the dis- 
crete ness of the particles accompanied by their associated gravitational polarization cloud (.Weinberell998l : lRostoker fc Rosenbluth 



1960|). This is shown by working out the Fokker-Planck equation for the evolution of actions of the particles in this random 



field. The potential of a mass m^ with action-angle variables J2, W2 on a particle 1 with action-angle variables Ji, wi is: 
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exp(i(ki-wi — k2-W2)) 



The fluctuating part of the potential created by the discreteness of the dressed particles is the sum over all particles 2 and all 
non- vanishing ki, k2 of potentials like (I39p . The rate of change of the action Ji of a particle 1 in this fluctuating fleld is: 

T Y^ Y^ Y^ ^ -1 exp(i(ki-wi -k2-W2) , , 

Ji = > > > Gmim2 Jki — — — — - — -— - • (40) 

^-^ ^-^ ^-^ I?kik2(Jl, J2,k2-rJ2) 

2 ki/0k27^0 

The braking and diffusion coefficients of the corresponding Fokker-Planck equation are obtained from, respectively, the first 
and second moments of the random change AJi suffered by the particle 1 in a time At. The averaging is performed on the 
values of the angle variables of particles 2 and on their action distribution functions. Equation (|38p is recovered that way. In 
th e calcu l ation of the braking coefficient, small departures from uniform angular motion should be accounted for, as shown 
bv lEckeiJ (1973) in ^ similar context. The Fokker-Planck form of equation (|38p . although equivalent to it, looks more complex 



than equation (|38|) itself because the braking coefficient involves the derivative of a Dirac distribution. 

5.3 Accounting for strong collisions 

Equation (|38|l and its quasi-homogeneous limit, equation (|47p . both result from a weak collision theory. Strong collisions 
involving substantial deviation of at least one of the colliding particles are not adequately described. This inappropriate 
description of the rare strong collisions can be fixed by limiting the range of impact parameters to values larger than some 
critical limit fe°^ which depends on the masses of the colliding species. This critical impact parameter for particles of species 
a and h is such that the typical kinetic energy in their relative motion be equal to their interaction energy, that is: 

GM maVHh _ GmaTrii, 

R rua + mb ^ bi^ 
Here M is the total system's mass and R a typical global size of it. Were this cut to be omitted, the expressions of the 
coefficients in equations (|38p and (|47p would diverge logarithmically at large wavenumbers, where the response function e 
approaches unity. This divergence results from the neglect of large deviations in strong collisions. A physically sound result is 
obtained by limiting the K integration in equation (|47|l to the domain | K | < K'^,. where: 



rab 



2-K 2-K M 



K^^ = ^ = :il . (42) 

b?,*! R rria+mt ^ ' 

Similarly the summations on the angle Fourier variables k^ (i = 1 or 2) in equation (|38|l should be limited, in the term 
associated to species a and 6, to values such that the physical wavenumbers along the quasi-intersecting orbits be smaller than 
ii'cr • This modulus of the physical wavenumber can be crudely related to the dimensionless angle wavenumber hy K — k/R, 
where 7? is a typical global size of the system and k the modulus of the angle Fourier variable. Thus, the summation on ki 
and k2 in equation (|38|l should be limited to wave vectors, the modulus of which is bounded by: 



,, , 2nM , , 

lk,l< . 43 

ma + nib 

When solving equation (|38p . the secular evolution of the response matrix e, the system's collective potential U{r, t) and 
the Fourier transform coefficients ?/>£ (J) should be followed in time together with the 1-body distributions. We return to this, 
in the case of spherical potentials, in section [S] Prior to that, let us discuss various limits and approximate forms of equation 
l|38[l and show that, as it should, it satisfies an H-theorem. The irreversibility stemms from the fact that information is lost 
when the real issues of collisions are replaced in the equation by average ones, in particular by averaging over the angles of 
the colliding particles. 



6 LIMITING CASES 

6.1 Homogeneous limit 

Although the limit of an homogeneous medium cannot be rigourously taken for a self gravitating system, it is nevertheless 
possible to assume local homogeneity at the price of artificially limiting the interaction distance between particles by cutting 
it at some characteristic size of the system. So doing, the effects of the collective dressing of the particles are still retained, 
albeit less precisely, but the effects of the structure of the system are only sketchily accounted for. 

In this limit, the system is regarded as homogeneously flUing a l arge cubix box of side L, on the surface of which 
periodic boundary conditions apply. This is the geometry considered bv I Weinberg ( 19931 ). Due to the assumed homogeneity. 
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the collective force F*^ vanishes and the unperturbed motion is rectilinear and uniform, whatever the state of relaxation of 
the system. The action variables are then proportional to the components of the momentum p and the angle variables are 
proportional to the components of the position r. Since the angles must be variables of period 27r, the angle vector must be 
w = 2tiy/L, which implies that the action vector is J = Lp/27r. The angle Fourier vector is k = LK/27r where K is the 
usual wave vector of Fourier transforms with respect to position. The frequency vector fJ is 2-kw /L, so that k ■ fi = K ■ v. 
The density-potential basis consists of functions proportional to complex exponentials, like exp(iK-r). A given element of 
the basis, a say, is characterized by its wave vector K. This can be accounted for in the notation by writing this wave vector 
as Kq, the corresponding angle wave vector being noted k^. The density function and the potential of the element a of the 
basis are both proportional to exp(JKc,-r). Their normalization factor must be such that the biorthogonality relation H22|) 
be satisfied, the density D°'{r) and the potential '(/''^(r) being related by equation H21I) . These constraints result in: 

^ ' 2V^LV2 ^ ^'^^ Li/2 |K,L| ^ ' 

The V'k's are the Fourier transforms of i/'"(r) with respect to the angles w, namely: 
,ti 2v^ (5(k„ - k) 

where S(ka: — k) is a triple Kronecker symbol. In this case the tp^'s do not depend on the actions and remain fixed while the 
relaxation proceeds. The response matrix e, calculated from its definition (|34p . is diagonal, its element aa being given, for iu 
in the upper half complex plane, by: 

' (^) = ^-2^ -JK^ JdP ,,K.., • (46) 

For real u, a -|-iO should be added to the singular denominator. Since a enters this relation by its wave vector Kc, e""{oj) can 
be regarded as a function e(Ka, c^j), or, more generally, as a function of a wave vector K and of a frequency uu. Because of the 
diagonality of the response matrix e and the simplicity of equation (|45p . the writing of the kinetic equation (|38p simplifies to: 



9tr(p) = E A'p' ^p- {Qa,ip,p') ■ (Vp-Vpor(p)/(p')) , 



(47) 



where the tensor Q^;, is defined by: 

^ ( M or^2 2 2 f,3r^ Sg ^(K-(V-VO) 

Qa6(P,P) = 2G m„mb dK — — ■ (48) 



J K* |e(K,K-v)|2 

Equations (|47p - (|48p are identical to equation (29) of lWeinberd ( 1993h when the quasi homogeneity of the system and associated 



absence of collective and external forces are accounted for. Equation (|47p can be written explicitly as a Fokker-Planck equation 
in the form: 

9t./"(p) = - Vp ■ (A„,r(p)) + i Vpv;: (s.r(p)) , (49) 

where the momentum drag and diffusion coefficients are: 

A.(P) = E AV/(P') ((Vp - VpO • ^„,(p, p')) , (50) 

b •' 

l.(p) = 2^ /dV/(p') ^„,(P,P') • (51) 



For electrical instead of gravitational interactions, the gravitational constant G should be replaced by l/47reo in MKSA units, 
to being the dielectric permittivity of vacuum. The electric force between like charges being repulsive instead of attractive, 
the minus sign before the second term of equation (|46|l should be changed to a positive sign and the masses replaced by 
the char ges of the particles. Equ ation (|47p then reduces to the Balescu-Lenard equation for homogeneous and multispecies 



plasmas (|Babuel Pevrissadll974l '). It implicitly accounts for the screening effect, which is embodied in the dielectric function. 
The integral on wave vector space in equation (148 P then need not be cut at small wave vectors because | £(K, K- v) | diverges 
as I K I approaches zero. For self gravitational systems the small | K | limit is unphysical, due to the absence of screening. 
The distance between interacting particles is limited in this case by the inhomogeneity of the system, a feature which is lost 
in the local approximation. If one were to insist on the quasi-homogeneous approximation, the integration over wave vectors 
in equation (|48p would have to be artificially limited from below to some minimum modulus Kmin ~ 2-k / R, where 7? is a 
characteristic size of the system. Little would then be gained over a more traditional Fokker-Planck approximation, but for 



the fact that equation (|47p still accounts for the collective dressing of the colliding particles. 
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6.2 Non-collective homogeneous limit 

When these collective effects are themselves neglected, which amounts to take £ = 1 in equation (|48|l , the usual local Fokker- 
Planck equation (|49|l is recovered, with braking and diffusion coefficients given by expressions (|50|l and (|5H) . e now supposedly 
being equal to unity. As above, the integral on wavevectors in equation (|48|l should limited to a lower cutoff at | K|= Kmin, 
to account for the finite size of the system, and to an upper cutoff | K j= Kmax, to account for strong collisions (section l5.3p . 
The coulomb logarithm is In A, where A — Kmax /Kmin- When e equals unity, the integration over wave vectors in equation 
H48[) can easily be performed. The result, which involves the relative velocity of the colliding particles g = v — v', is: 

A,(p) = -47rG' ln(A) ^ mamb{ma + mt) /rfV/(p') 4 ' (^2) 

l„(p) = +47vG^ ln(A) Y, ^Wb fdVf'ip) 1^1^ , (53) 



where I is the unit second rank tensor. This is identical to the Fokker-Planck equation presented, for example, in lBinnev fc Tremaine 
(1987)), equations (8A.10). The collective dressing becomes important when | e |~^ in equation (|48|l largely differs from unity. 



As shown bv lWeinberg ([1993), this happens when the system is not far from being unstable, for example when its size becomes 
of order of the Jeans length, the complex zeroes of £(K,aj) lying close, but below, the real axis. 



7 AN H THEOREM 

Equation (|38|) satisfies an H theorem which states that the statistical entropy: 



s(t) = -^ /d'wi/"d'jir(wi,Ji,t)in(r(wi,Ji,f)) 



(54) 



increases with time. Because the relaxing distribution functions depend on actions only, the integral over angles reduces to a 
mere multiplication by a factor 87r"^, so that: 

0^(1) _ ^^3 \ ' ; jsr / 1 I i„ rf^ri. +^^ \ / ^ ^"i't_ +\ \ ('55') 



dt 



-8^'^ foi (l + ln(r(Ji,t))) {dtr{Ji,t) 



The time derivative of /" is given by equation ()38p which can be symmetrized by substituting to the first operator ki-Vj^ the 
operator ki-Vj^ — k2-Vj2. The contribution associated with the added operator k2-Vj2 vanishes on integration over J2. This 
can be seen by using the flux divergence theorem in action space, recognizing that the surface integral over the boundary of 
the physical J2 domain vanishes. Indeed, the expression on the right of the first operator ki- Vj^ in equation ((38} represents, 
up to its sign, the flux in action space at Ji caused by collisions with particles having action J2 or the flux at J2 caused by 
collisions with particles of action Ji. The physical domain is limited in action space by a boundary at a finite distance and 
extends to infinity in certain directions. The flux through the boundary at finite distance vanishes because the action vector 
of no particle can evolve through this boundary from the physical to the unphysical domain. The fiux at infinity vanishes 
because /''(J2) decreases fast enough. This justifies the above-suggested substitution. The expression of 9t/"(l) given by 
equation (|38|) . modified as described, when inserted in equation (|55p . gives the following expression for ds/dt: 

ds{t) 






dt 

••• (ki-Vji-k2-Vj2)<5(ki-rii -k2-r22)|Okik2(Ji,J2,ki-ni)l~'(ki-Vji-k2-Vjjr(Ji)/(J2) ■ (56) 

This expression is further symmetrized by combining it with the equivalent expression obtained by exchanging species indices 
a and b, actions Ji and J2 and Fourier variables ki and k2. The resulting expression is then integrated by parts, using the 
fiux divergence theorem in either Ji or J2 space. As explained above, the contribution of the integral on the boundary of the 
action domain or at infinity vanishes. We are eventually left with the positive expression: 

^ - +32.^ E E E E h^ l^'J^ ^'-i-i ^(k..».-k2.»2) ((k.v.. k2.v..)r(Jo/^j2))- 

'^' a,r^i:J ' |Ok,k.(Jl,J2,kl.fll)l' ,f"(Jl)/^(J2) 

This establishes that the statistical entropy of the system is monoto nically increasing. Since th e entropy of a self gravitating 



system of a given total mass and energy is not bounded from above (JBinnev fc Tremaindl 19871 ). the increase of the statistical 
entropy does not lead, as in homogeneous gases or plasmas, to a state of thermodynamic equilibrium. When the system become s 
sufficiently centrally condensed, a gravothermal instability develops ( HenorJll96ll : lAntonovlll962l : lLvnden-Bell fc Woodlll968l ). 
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The quenching of this instability by the formation of binaries is not described by equation (|38p . because the formation of 
binary systems results from triple collisions (that is, from third order correlations) and is of the strong interaction type. 



8 EVOLUTION OF A SPHERICAL POTENTIAL AND BASIS FOURIER COEFFICIENTS 

When the relaxation proceeds, the distribution functions /"(J, t) evolve according to equations (|38p . This causes a slow secular 
change in the average potential U{v,t) and in the response matrix e{ijj) (equation (I34[l 'l. The basis potential functions '4>°'{r) 
are time-independent, but their Fourier transforms with respect to angles w are not because they depend on the orbits of 
the particles, which slowly change with the potential. The Fourier coefficients V'k(J) S'l'e then indirectly related to the slowly 
evolving potential U{v,t). 

Thus, equation (|38[) is not an autonomous equation for the distribution functions /"^(Jjt). The response matrix e{uj) of 
equation (I34|l is a functional of those, which also depends on the angle Fourier transforms il}^{J,t) as does the quantity V 
present in equations (|35|l and (I38p . The kinetic equation (|38p must then be completed with equations describing the evolution 
in time of the average potential U{r,t) an d of the a n gle F ourier coefficients ip^{J,t) of the basis potentials. This aspect of 
the system's evolution is not considered bv IChavanisI (2007)). In this section, the time t will be restaured, though only where 



necessary, in the list of arguments of functions. The potential U{r, t) derives from the mass density: 

D{v, t) = Yma fd'p r (r, p, t) • (58) 



The corresponding gravitational potential is obtained from its expansion on the biorthogonal density-potential basis by 
equation (|24|l . From appendix |B| it is found that its partial time-derivative is: 

dtU{r,t)=-J2»^''Gmai'°'{r)fd'j 9t(r(J,t) (V'o(J,t))*) • (59) 

a 

Equation l|59p . as well as equations (I34|l . pSfl and (|38|l . call for an equation for the time-evolution of the coefficients tp^{J,t): 
V'£(J,i)= /d3we-*-"'i/;"(r(w,J)) ■ (60) 



An explicit expression for these coefficients when the potential is spherical is derived in appendix [B] which also gives a 
summary of angle and action variables for a particle moving in a spherical potential. In this case, the coefficients tp^{J,t) 
vanish when the wave vector k has non-vanishing ^2 or fca components. The non-vanishing coefficients depend only on the 
radial fci component, hereafter noted k. For conciseness, the variables J, which are mere parameters, are omitted wherever 
possible. After some calculations, we find that: 

:c,, „ 2^ , ^ r^''^ COS Wk(r, t) iP°'(r)dr , t,w x T kUAt) dr' 



.(t) y/2{E{t)^Uir,t))~Ji/r^ Jr.it) y/2{E(t) ^ U{r',t)) ~ Jl/r 

When the potential changes, the radial distances rp and va of the periapse and apoapse change accordingly: the bounds of the 
integrals in equations (|6ip . (|B5|) and (jB2|) are time-dependent. These integrals are singular, though convergent, the periapse 
and apoapse distances being the zeroes of the square root denominator: 

q{r,t) = 2{Eit)-U(r,t))-j'i/r^ ■ (62) 

These zeroes are simple when the orbit is not circular and merge into a double zero when it is. This latter situation can be 
dealt with by a limit process, in which simple zeroes rp and va are made to converge to eachother. We then assume that rp 
and Va are simple zeroes. An index P or A denotes the value of a function of r at rp or rA respectively. Integrals like: 

Jrp(t) \/q{r,t) 

or similar ones can be expressed in terms of a variable ^, the values of which remain constant at the changing periapse and 
apoapse. This variable is defined by: 

r = rp{t) + ^{rA{t)-rp{t)) ■ (64) 

To each of these two types of radial variables, r or 5, a time variable, t or r, can be associated, it being meant that t = t. This 
introduces two types of time derivatives: dt, which is at constant r, and Ot, which is at constant ^. Ordinary time derivatives 
with respect to r and t are identical and are denoted by a dot. Partial time derivatives with respect to r and t differ and are 
related by: 
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dr^dt + 



(rJ^-I^)+r,(l^^^))dr ■ (65) 

V \rA~rpJ Xta — rp J J 

The partial derivatives with respect to r and ^ are simply proportional: d^ — {ta — rp) dr- At the periapse or apoapse the 
function q{r,t) vanishes, whatever r. Hence, drQ = at these points. Differentiating the equation q(r,t) — 0, rp and va are 
found: 

dtq(rp,t) . dtq[rA,t) 

^P^- ^— 7 TT rA = - j—f -r ■ (66) 

drq[rp,t) drq(rA,t) 

The partial derivatives of q(r,t) (equation (|62[) ') are: 

dtq{T,t) = 2(E~dtU{T,t)) , drqir,t) = 2(^^drUir,t)j ■ (67) 

It can be checked from equations (|66|) and (|65p that [dTq){r-p,t) — dTq{rA,t) — 0. Equation H65|) implies that: 

d^r = rp(^:^^^)+rA(^^^^) , a.g = 9.g + (fp (^lA^) + f^^-I^) ) 9.g • (68) 

\rA — rpJ \rA — rpJ \ \rA—rpJ \rA — rpJ J 

The time-derivative of I(t) (equation H63p ) is found by changing the variable r to ^: 
^ ^ / rA-rp \ ^ ^ r^W m(r,t)dr / 9.m _ 1 ^^\ 

It is important to note that the last term in the parenthesis of the integral in equation (|69p is regular since the numerator, 
drq, vanishes at rp and ta, where q{r,t) does. The right hand side of equation (|69|l then consists of convergent integrals. 
When this method is used to calculate the time derivatives of E{t) and ^i{t) from equations (|B2|) and (|B5|) . the following 
results are obtained: 

E ^ ^ r' ^r ^m , n. ^ -n. (tA^) + ^ r' -^ ^ • (70) 



-^ '- - >/g(r) ' \rA~rp) ^n J^^^^^ Jq{r,t) ?('■-*) 



'rp(t) 

The same method is used to calculate drWk'- 

"i \rA-rpJ 2 Jrp(t)^/^iy7^ 9('^'*) 

The angle Fourier coefficient i/'^(J,i) is given by equation (|61|l . which is of a form similar to equation ([63}. Using the general 
result (|69|l . the time derivative of V'fc (Jji) is found to be: 

^nJ,t) = gi ^," + (tA^ ^^.sn^n. M^'r^ (su.W,ir,t) a.W, - cosl^.(.,t) ("t^ _ ^)) . (72) 

Equation (|72p describes the time-evolution of the Fourier coefficients i/jfc(J,i). The auxiliary r-derivatives which enter this 
equation are given by equations (|71l) . (17011 . (I68p . (|66p and (|67l) . The other quantities entering equation (|72p are expressed in 
terms of the potential by equations (|61l) . (|62l) . (|B4p and (IBSp . All these relations eventually link the time-evolution of f/;^ to 
the time-evolution of U{r,t), which is itself described by equation (1591) . 



Equations ((38]), (|59p and (|72p form the system of coupled equations for the distribution functions /"(J,i), the average 
potential U{r,t) and the angle Fourier coefficients tp^{J,t) that we have been seeking for in this section. This system involves 
the auxiliary equations mentioned above, as well as equations (p4)) - (p5]) . 



9 CONCLUSIONS 

Kinetic equations for the coUisional evolution of the constituents of self-gravitating inhomogeneous systems have been derived. 
These equations (|38p surpass in consistency the usual Fokker-Planck equations (|49p - (|52p - (|53p . The latter are unsatisfactory 
from a principle point of view, being local and non-collective. 

By contrast, the proposed equations fully account for the system's inhomogeneity and for the collective gravitational 
dressing of the colliding particles. 



Equations (|38p describe the evolution of distribution functions in action and angle space, which is possible when the 
hamiltonian associated with the average potential is integrable. 

Physically, these equations describe the evolution of the distribution functions in action space as a result of the weak 
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gravitational noise caused by the discreteness of the particles, dressed with the polarization clouds that their own gravity 
induces in the system. This gravitational polarization is accounted for in equation (|38p in a manner that is fully consistent 
with the distribution functions, as they are at the moment. 



9.1 Properties of the kinetic equation 

Equation (|38p is the sum of a second order derivative term with respect to actions and of a first order one. It therefore basically 
is of the Fokker-Planck type, although it is definitely simpler in the form of expression (|38p . The diffusion coefficient involved 
depends on the 1-body distributions themselves, in particular through the factor \T}\^^ which represents the effect of the 
dressing of the colliding particles by the gravitational polarization induced around them by their own influence. 

Unlike in electrical plasmas, the polarization dressing in self-gravitational systems does not cause any screening of the 
interaction, which remains effective even between distant particles. The mutual distance of such particles is limited only by 
the finite size of the system. Were the gravitational influence of particles on their surrounding to be neglected, the response 
matrix e (equation (|34|) ) would reduce to unity and the coefficients of the corresponding Fokker-Planck kinetic equation would 
simply be averages by the distribution functions of functions of velocity, as in equations (|52[) - (|53|) . 



It is apparent from the developments of appendix [X] which lead to equation (|38|) . that the k component in angle Fourier 
space of the gravitational polarization response given to a particle has frequency u — k-fJ. This means that the polarization 
cloud which accompanies a particle forms a structure in angle space which vary as w — i7i: it corotates in angle with that 
particle. 



The presence of the Dirac function 5(ki ■ fii — k2 • ^2) in equation (|38p indicates that particles interact resonantly. 
This certainly is an important physical property of remote interactions, for which the components of the angle wave vectors 
ki and k2 must be small. For closer encounters, the modulus of these wave vectors is larger and the resonance condition 
ki • fJi = k2 • fl2 becomes less selective, being more easily satisfied. 

The correlation function has been calculated on the basis of a linearized theory, which is justified by the weakness of 
the average interactions in this many-body system. This means that the trajectories of the particles during the collision are 
regarded as being the unperturbed trajectories. Similarly, the gravitational polarization cloud around any one of the colliding 
particles is calculated as if the partner in the collision were not present: equation (|38p is still a weak collision approximation. 
A cutoff at small impact parameters is therefore needed to account for the rare strong collisions. 



Equation p8|) takes full account of the inhomogeneity of the system, which is embodied in the dependence of the distri- 
bution functions on the actions J's. It requires no artificial cutoff at large impact parameters. The details of the trajectories 
followed by the particles in the present gravitational potential are also fully accounted for, being implicit in the relations 
which link the angle and action variables to the position and momentum ones. These relations depend on the actual global 
gravitational potential of the system, which slowly evolves in time together with the distribution functions. 

The density-potential basis functions il)" (r) are choosen at the beginning of the calculation once and for all, but their angle 
Fourier transforms i/;£(J), which depend on the actual trajectories of the particles, change with time because the trajectory 
of a particle of given actions slowly evolves with the general potential of the system as the relaxation proceeds. As long as it 
suffers no collision, a given particle keeps its vector J fixed because the actions are adiabatic invariants. Collisions, however, 
cause a secular evolution of the functions /"(J), which is exactly what equation (|38[) describes. 

The description of particle motions is made simple by the use of action and angle variables. Their complexity is embodied 
in the supposedly known relation between position and momentum variables and action and angle variables. The usefulness 
of equation I 38ll is therefore limited to systems for which this relation can be established, either analytically or, possibly, 
numerically 



Pichon fc Cannonlll997l : iMcMillan fc Binnevlliooi ') 



While the relaxation proceeds, the gravitational potential and the orbits of the particles evolve. As a result, the kinetic 
equation must be completed by evolution equations for the potential and for other relevant quantities. Section [5] establishes, 
for spherical systems, this set of coupled equations. 
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APPENDIX A: DERIVATION OF THE KINETIC EQUATION 



Equation (|37(l must be expressed in terms of angle and action variables, using the adopted density-potential basis. The 
integration on the dynamical state of the particle of species h should be carried out first. The integral over the variables 2 in 
the second line of eQ. (|37|) is similar to equation (|15l) . with A — Gniamb, M{2) — C/^(2, 2',aj') and N{1) — 5p(l, I'jti^) and is 
thus expressed in angle and action variables by equation (I27p . Since M(2) is as in equation (|13p . the a^ coefficients are those 
of equation (|28|l . the species indices being now h and q instead of c and p and the parameters being 2',aj' instead of 1' ,ui. 
This leads to the following change in equation (|37|) : 



^ /'d2 (F<,,(1, 2) • Vp, ) ^^"(1, 1', c^) g,*(2, 2', J) 



i Grngmg y^ y^ (e ^(t^'))° i>t*{J'2 



pi(k2-wi-k2-W2) 



k2 k' 



((^£,(Ji)*k2-Vj, -(Vj,V£2(Ji))-Vwi) g;(l,l',a-)) ■ (Al) 



Note that, as a general rule, operators act on everything on their right, up to the end of the expression or to a closing delimiter. 
Using the relation (|A1|I . the collision operator (|37|l can be written as: 



C"(/) 



E 



dij' 



dr dl' d2' n^ / ^e-""+">5:5: ^Gm„m,(s-(a.')rVfr(J'.; 



, 2n 



k2 ki 



„i(k2'Wi-k2-W2) 
CJ' - k^ ■ fl'^ 



^£,(Ji) ik2 • Vj, - (Vj,V£2(Ji)) ■ Vw,)e,"(l, l',c.) ) (f,,(1',2') . (Vp, - Vp,) r(l')r(2')) 



(A2) 



The 1-body propagator ^/^(l, l',a;) is then Fourier-expanded with respect to both angles wi and wi according to equation 
(|19[) and this expansion is inserted in equation (|A2p . It then appears that C°'{f) depends on wi as exp(i(ki + k2) • wi). Since 
during relaxation /"(I) remains a function of Ji only, it is possible to average over wi without loss of information, which 
brings a Kronecker factor 5(ki + k2), such that ki — — k2 = k. The angle-averaged form of the collision operator is: 



cw - -EEf-Z-'Z-'/JI^ 






r.^(\){e-\J)) 



.'\\al3 



\J Gniar 



<:(2' 



^ • nr '"''"' E^'"'''™^ Grk;(i.i' ^) F,,(i',2')-(Vp, -Vp;)r(i')r(2') 



(A3) 



The following calculations are somewhat similar to those carried out for an homogeneous plasma by llchimarul ( 19731 ). The 
expression (|A3|l can be split into two parts, one, Ci{f), being associated with the operator Vp' in the last parenthesis and 
the other, C2{f), being associated with the operator Vp' , so that: 

CV)-Cr{f)+C,V) ■ (A4) 

The terms Fp,(l',2') -Vp' .r{2') and Fp,(l',2') -Vp/ /''(I'), multiplied by other functions of 1' and 2' respectively, are subject 
to an integration over these variables. The structure of these expressions being similar to equation (|15p . they are expressed 
as in equation (|18p . noting that f^{l') and ,f''{2') do not depend on angles. The coefficients aa of the development on the 
density-potential basis (equation (|24|) ') which enter equation (|18p are calculated from equation (|26|) . with appropriate M 
functions. Integration over angles w^ or W2 can then easily be carried out. We are left with: 



Ciif) = -i 



' EE P^i^l, ^'~^ '^''^''^ Y.^^^'Gfn..m,ml (k. v., ) 



r~^{i){e-\j)r^ 



?/''- '= (-,"")?/« 



r^t:{2'WL*{2')P{2') 



k2 ■ Si2 



(A5) 
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C2V) = +i 



L^Z^ 



dr 



— / — e ^+ 2^(87r G) mampm, (k-VjJ 

B Jb' ^. 



^-ui)ie-\^')r 



\^ U^T' n'^p 



d'j[ G"^,(i,i',^)^r(i')r(i') 



k' 



^'-K- ^'2 



(A6) 



The integrations over r and a;' which appear in equations (|A5|) - (|A6P are of the general form 

The integration over r is regular, and straightforward, when u + u}' has a negative imaginary part. Otherwise the result must 
be obtained by analytical continuation. This means that, whatever ui: 



(A7) 



K.)^l t^^/H.(-') 



(A8) 



The contour B' passes above all singularities of gij-o')- For a stable system these are all below or on the real axis. When uj is 
low enough in the lower half complex plane C~ for —a; to be above B' , the integration on uj' can be carried out by closing 
the contour B' at infinity in the upper half complex plane C^ , using the theorem of residues at the unique singularity in 
the closed up contour, which is at uj' — —u). The closing of B' in C^ is possible because in the present case (see equations 
(|A5I - IX6)) ) g{Lo')/{uj + u]') decreases at infinity as | uj' |~^, which means that for such oj's, h{ijj) = f{uj)g{—uj). Analytical 
continuation extends this result to other uj's. The two parts of the collision operator then reduce to: 



Ciif) =~^ I J^E (STv'Gfm^k-Vj,) 



^-Ui){e-\-o.)r^ 



J2"^^J2 



Jd'ji G^-,,(l,l',a;)V\,(l') (k;.Vj;r(l')) j I E "^' E /^'-^^ 



<;(2')V!;'(2')/^(2') 



(A9) 



l^_,.,^^«'9 



C.Kf)^+^l ^Y.^87v''Gfma{k-Vj,) V-k(l)(e-'(-c.)) 

Y^m^Y^hi G-,(i,i',a;)v.-(i')r(i')| ( E "^? E A''^^ v^,(2')<;(2') ^'I'^^'t^/^,^^ 



q K 



(AlO) 



The last parenthesis in the second line of equation (|A10|) is {S^^ -— e ' {~lo))/{9,-k''^G). Similarly, the last parenthesis in eq. 
dm is ii•''^(-LJ)/87r^ where the matrix //""(u) is defined by: 



Ji-"''H = 8^^E"^'E A'-^ 



, c;(j')v!k'(J')/^(J') 



- - k' • n(J') 



(All) 



q k' 

The double Fourier transform C^, of the propagator which is present in the first parentheses of equations (|A9|I - (|A10|I can 
be read from equation ((36]): 

gk-K(M'.^)= (Z^^f^T) (8^^p^('-+k;).(J.-j;)+ """^^l^^^^X)^'^^ ^^(Jr)(.-(..))N^-,(J'0) . (A12) 
Using this, the first parenthesis of the second line of equation (|A9|) . VJi(l,a;), can be written as: 

^n M \ * "^o (k-Vji./"'(l)) ,A/,x/ -1/ N\^7 /K-to\ 

87r'' IX) — k • iJi ^ ' 

The first parenthesis of the second line of equation (jAlOll . V^2{^i^)^ i^ similarly calculated and expressed in terms of the 
matrix H defined by equation (|A11|I : 



Inserting equation (|A13|I in equation HA9p we get: 



(k.Vj,r(l)) ^,^^^ (e-i(^))^"i/M7(+^) 



a; — k • fJi 



(A14) 
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(A15) 



r-.{l) {e-\-u.)y' V^^(l) {e~\+<.)Y^'H^-{-u.) ^^^^^ 



Inserting equation (jA14p in equation (|A10[) we get: 

( ^,r(i)'/^rji) + (k-vj,,r(i)) ^i{i){e-\+u^)r^H''-'i+o.) ) 

Gathering equations (|A15|) and (|A16|) . the following expression is obtained for the collision operator 



(A16) 



C"(/) 



/£EG'-^(k.V.J 



rk(l) 



Lj — k ■ ill 



+ {(e-\-u;)r-S"^)r-lW^rW - (k-Vj,r(l)) ^ptW{e-\+u;)f^H^-{+o.) 



'G 



(A17) 



The last term on the third line vanishes on integration over uj. Indeed, the Bromwich contour must pass over all singularities 
of the function /(oj) in eq. (|A7[) . that is, in equation (|A17|) . over all singularities of functions of +tj. The contour B can be 
closed at infinity in the upper complex plane, which gives, for the fourth term of the square bracket, a vanishing result. The 
two terms in the second line of equation (|A17|) can be associated, yielding the following expression for C"(/): 



C'U) 



ki 



^- 



.(1) ^^33^^ {{e-\-.)r-nr_i.i^)^-^) 



I'^T. G^-^(ki-VaJ r^.^il) ^^^1-^^nm is'\~c.)rr^^il)ie-\+o.)r[H^^-.)+H-^^+..)) 



(A18) 



To evaluate the second line of equation (jA18[l . the integration contour B may be lowered to the real axis. Rigourously, uj 
pertains the upper complex half plane and can only be consider real in a limit sense when the contour B descends to the real 
axis. Real singularities at a; = k ■ $7 must therefore be avoided by the contour by skirting them from above, so that: 
lip 



iTi5{uj — k • n) 



(A19) 



cj — k-Jl Lj — k-n + iO Lj — k-J7 

where V is the principle value distribution. Conversely, when a; descends to the real axis, —a; rises to it from below, so that: 
1 (-1) V 



-cj + k-n Lj-k-n-iO t^-k-n ^ ' 

The sum H'^'<{-u]) + H'<^{+uj) calculated in this limit is: 






ml rJii'wL (1') s{^ - k'l • fi'or (1') 



(A20) 



(A21) 



Thanks to the Dirac function in equation (|A2ip . the second line of equation (|A18|) is easily integrated over lj. Where conciseness 
demands, we note: 

cJi — ki'ili , cji — li-ifli ■ 

The first line of equation (|A18[I can be disposed of by closing the ui integration contour in the lower half complex plane, which 
is possible because the integrand declines fast enough at infinity. The system being supposedly stable, all the singularities of 
(£~^(— a;))™'' are in the upper half plane. The contour then encloses only the real singularity at a; = ki ■ rii and its sense 
brings a factor — 2i7r when using the residue theorem. The expression of the collision operator then becomes: 



C''{f) = iGml Y, (ki-VjJ 

ki 



rk,(i) {(e-\-u^i)r^ - r") ^!;^(i)/.(i) 



+ 



Y^Y. h'J'- ^-'^< {r-kA^){e-\-.[)Y%%,y)) (<(i)(e-i(+c.i))^^v^r(i') 



>,\^'^.,.i*,.'<\ krVj,(r(l),r(l')) 



cji — cji + iO 



(A22) 
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This expression is then symmetrized. Half of the term on the second line of eq. (|A22|) is added to half of the same expression, 
modified by changing ki into — ki and kj into — kj. This leaves it invariant, except for the last denominator, which is changed 
into — (oji — cji — iO). The half difference brings a contribution —inS{ijj[ — uji). The term on the first line of equation (|A22|) 
may be similarly symmetrized. When changing ki to — ki , the argument of the inverse response function changes sign. The 
change of the response function when the sign of its real frequency argument, cjr say, is modified may be found by noting that 
its real and imaginary parts, e' "^ and e" " , are: 



/Q/3/ N r 

e (UJr) = 



£""''K) = ^^ /"dVi 87r*Gm^ V'£;(1)<(1) (ki-Vji,r(l))5K-ki-ni) 



(A23) 

(A24) 



q kl 

The conjugation relations (12511 can be used to show that; 

e-^i-LOr) = (e-^i+UJr))* , ie-\~LOr)r^ = {{e-' {+OJr)f^ )* , (A25) 

where the basis element a associated with a is defined by equation (|23l) . Using this, the expression T, defined by: 



T = ^(ki.VjJ 



kl 



r-^,{i) {{e~\^\,rn^)r^-s'^^) ^!yi) ,r(i) 



which is present in the first line of equation (|A22|) is symmetrized to: 
2 



T=-Y. |(ki-Vai)^£,(l) ({s-'{u^)f'-({e'\ur)) 



<(i)r(i) 



(A26) 



(A27) 



kl 



As expected ( Nelson fc Tremainelll999l ). this expression involves the antihermitian part of the matrix e ^, which may be 
expressed in terms of the antihermitian part of the matrix e as: 

{s-')-{s-y=s-\s^^s){er' ■ (A28) 

The matrix e^ — e is calculated from equations (|A23I) - (IA24|) : 



(e''"H)*-£"''H = -^^^ WTv^Gml fd'ji 5iLo--k[-n[) (k'l-Vj- J'(f')) ^^(0 <; (l') 

The term T in equation (|A27I) can then be written as: 

r = 8J7r*^^^ d^j[ Gml 5(ki- fli -k'l- n'l) (kl- Vji 

n W, 1,' ■' 



q kl k; 



<(l)(e-^(c.O)'''<;(l') 



(k;.vvr(i')) 



(A29) 



(A30) 



and inserted in the first line of equation (|A22I) . The square modulus factor in equation (|A30|) is | V^^^^^i {Ji,J'i,uji) \~^ 
(equation ((35])). The second line of equation (|A22|) can be treated similarly. From equation p5[) . one of the parentheses is 
(^kik' (Jij Ji)'-^i))~^ s-'^d the other is its complex conjugate, which can be shown by using the conjugation relation (|25p . 
When all these symmetrizations and substitutions are made, the collision operator is finally written as: 



q kl k' 



C"(/) = X1EE/ "^'^1 STT^G'mX ki-Vji 



5(ki-ni -k'l-fi'i) 



2?kik;(Ji,Ji,ki-ni; 



ki.vji-k'i-Vj; r(Ji)r(j;) 



(A31) 



APPENDIX B: VARIABLES AND FOURIER COEFFICIENTS FOR SPHERICAL POTENTIALS 

Bl Angle and action variables for a spherically symmetric potential 

The motion of a particle in a spherically symmetric potential is best described in spherical coordinates r, 6, tp, the variable 
r being the distance to the center, 9 the colatitude measured from the pole associated with the coordinate polar axis z and 
(p the azimuth measured from an arbitrarily defined origin in the equatorial plane. Let U{r) be the gravitational potential, 
an increasing function of r approaching at infinity, which is provisionnally treated as constant in time. The fact that U{r) 
actually slowly evolves as the relaxation proceeds is addressed in section |8] Without loss of generality, the particle may be 
assumed to be of unit mass. A dot indicating time derivative, the conjugate momenta to r, 6, ip are: 



Pr 



pe 



p^ — r sin 6 ip 



(Bl) 
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Figure Bl. Angular parameters associated with the projection M of the particle on the unit sphere. O is the origin of the spherical 
coordinates, N the origin of the colatitudes and V the ascending node of the orbit. The orbital plane is OVM. 



In a constant potential, the energy E oi a. particle is a first integral, as is the vectorial angular momentum L, that is, its 
modulus L, its projection Lz on the polar axis and the direction of its projection onto the equatorial plane. The angle and 
action variables are deduced from the variables r, 6, ip, Pr, Pe, Pu> by a canonical transformation, the generating function of 
which is the solution to the Hamilton- Jacobi equation. iGoldsteinI ( 19561 ) shows how to construct angle and action variables 
Wi, W2, W3, Ji, J2, Ji in the case of a newtonian potential. Simil ar results for a general spherical potential are also well known. 
They can be found, for example, in Tremaine fc Weinberg (1984)1) or in Saha (1991). A summary is presented in this appendix. 

The orbit in a spherically symmetric potential being plane, the periods of the azimuthal and latitudinal motions are equal. 
This introduces some freedom in defining the actions, which can be taken advantage of to impose that one of the angles, W3 
say, be a first integral, associated with the direction of the equatorial projection of the angular momentum. The origin of the 
constant angle ws can be chosen to coincide with the origin of the azimuths and the origin of the radial angle variable w\ 
may be placed at some fiducial periapse. The angle and action variables then are given by the following expressions: 

y/2{E - U{t)) - J|/r2 dr , J2 ^ L = {pi + p%/ sin^ ef'\ 

I rp 

W2 = 1p ^ 



Jl = - 



Wl 



J'6 



■PV' 



(B2) 



w^ = (fi — arcsin (cot 6 cot /?) ■ (B3) 



y/2{E - U(r')) - J!/r'^ Jp y/2(E - U{r')) - Jl/r'^ 

P represents the position of a particle passing at the point M — r,9, ip, with momenta pr, pe, Pip when it reaches a fiducial 
periapse of its orbit. For given position and momenta, the action and angle variables in equations (|B2|) - (|B3|) depend on the 
potential U{r). The radii rp and ta are the distances to the origin of the periapses and the apoapses of the orbit of a particle 
with actions J = (Ji, J2, J3). They are given by the equation: 



2(£-t/(r))-j|/r^ 







(B4) 



and they depend on E and J2, that is, on Ji and J2 but not on J3. There being many different periapses, the fiducial one must 
be defined not only by its spatial location, but also by the time at which the particle passes there. The sign ± in equations 
(|B3|I should be taken as + when the particle visits the fiducial periapse P before it passes at M and — otherwise, fii and 
Q2 are the pulsations of the radial and latitudinal motions respectively (equations (IB5p ). The angle ^ is the azimuth of the 
present particle's position in the orbital plane, measured from the ascending node (figure iBlj) . The angle W2, which varies 
linearly in time, is the mean angular motion of the particle in the plane of its orbit. The constant angle UI3 is the azimuth of 
the ascending node in the equatorial plane. The angles wi and W2 are expressed as radial integrals following the sense of the 
particle's motion, whence the presence of an absolute value of the differential element in equations (|B3|) . The boundaries of 
these integrals on r' have not been written as rp and r because, depending on the relative position of the particle and the 
fiducial periapse, the integral may be extended over several successive senses of the radial oscillations. The ratio J3/J2 is the 
cosine of an inclination angle /3 (figure IbT|) defined by cos/3 — J3/J2. The latitude of the particle oscillates between ±/3. The 
frequency ils vanishes and the frequencies ili and 0,2 are given by: 



TV 

Til 



dr 



02 



Oi 

TV 



.h 



dr 



(B5) 



^2{E - U{r)) - Ji/r^ ' ^ J^^ '"' ^2{E - U{r)) - Ji/r^ 

They are both positive. The angle variables mi and W2 then increase linearly with time with the frequencies Oi and ^2, 
changing by 27r in a complete, respectively radial and latitudinal, oscillation. Equations (|B2[| - (|B3|) give the angle and action 
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variables in terms of the position and momentum variables. These relations may be inverted to give the latter in terms of the 
former. This however involves the inversion of the implicit relation (|B2[) to obtain _E as a function of Ji and J2 and of the 
first of equations (|B3I) to obtain r as a function of wi, Ji, J2. 



B2 Basis Fourier coefficients for spherical potentials 

The basis expansion coefficients which correspond to the density distribution (|58|) are obtained from equations (I17p and 
When, as here, the distribution functions do not depend on angles, their integrals over angles in equation (|26p are simply 
proportional to the k = Fourier coefficient of ^p"" , or equivalently of ^p°' (equation (|23p '). which is the complex conjugate of 
"00 (J)i) (equation (|25l) ). This coefficient depends on time, due to the slow variation of the orbits. We find that: 

ac{t) = -^ Sn'^ma d'j /"(J,*) (V'o(J,i))* • (B6) 

a 

From equations (|2ip and (|24l) . the gravitational potential is U{v, t) = G aa{t) tp"{r), a being a dummy index. Its partial time 
derivative is given by equation (|59|l . We also need some explicit expression for the angle Fourier coefficients tp^{J,t) We also 
need some explicit expression for the angle Fourier coefficients ^^{J,t) of the basis potentials (equation (|60|) V The relation 
of the position r to the angle and action variables w and J depends on the potential U{r, t), and thus on t. One could think 
of evaluating i/'£(J,i) for a given potential U{v,t) by just calculating the integral over angles in equation (|60[1 . The position 
vector r would then have to be expressed in terms of the angle vector w, for given actions. This cannot be done explicitly in 
general, since the relations (|B2[l - (jB3|) would have to be inverted. For spherical potentials, it is easier to change the variables 
of integration wi, W2, ws in equation (|60|) for position-type variables r, i/), ws (figure iBlj) . For given actions J, these variables 
are related by the equations (IB3|) which can also be written, with the notations explained above, as wi = Wi{J,M{r),t) and 
«;2 = V + W^2(J,M(r),t), where: 

Jp ^2{E{t)-U{r',t))-Ji/r'^ Jp ^2{E(i) ~U{r' ,t)) - J?Jr-^ 

The jacobian of the transformation from wi, W2, W3 to r, tp, 103 is | dWi/dr \. For a spherical potential, the basis potential 
functions can be chosen to depend only on the radial distance r. Equation (I60p then becomes: 

t/<£(J,f)= /|dr| / dv/ dws ^ "iWV' (r) ^ exp{^i{kiWiiM{r),t) + k2W2iM{r),t) + k2iP + ksWs)) ■ (B8) 



/o y/2{E-U{r,t))-Ji/r^ 

The integrations over the angles Tp and ws reduce to Atv'^SI 5? where the S's are Kronecker symbols. Thus the coefficients 
i/)£ differ from zero only when the ^2 and fca components vanish. The radial integration is over a complete oscillatory cycle of 
the variable r, from rp to rA and back. The coefficients ■i/)^(J,t) depend on the potential U{r,t) and on the fci component of 
k, hereafter simply noted k. Equation (|B8I) then reduces to: 



V'fc(J) = 47r^ni(t) i'ldrl —i-- ^' • (B9) 



The cycle integral over r in equation (|B9I) can be separated into an ascending part, in which r increases from rp to ta, and 
a descending part in which it decreases from rA to rp. Let W^{r) be the value of Wi{M{r),t) during the ascending part 
and VKj~(r) its value during the descending part. Wi{M{r),t) is a monotonically increasing function along the oscillation. Its 
value at the apoapse is tv. Equation (|B7|) shows that tv — W^{r) — W^{r) — tv. Defining Wk{r,t) — kW^{r,t), equation (|B9|) 
is turned into equation (|6ip . 
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